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1.  Introduction 


Kanel  et  al.  (1992)  shock  loaded  K-19  glass  in  a  normal  plate-on-plate  impact  test,  and  the 
VISAR  (velocity  interferometer  system  for  any  reflector)  measurement  of  normal  velocity  at  the 
free  surface  contained  a  second  plateau  that  they  interpreted  as  evidence  of  a  failure  wave 
(figure  la).  Impact  by  the  flyer  plate  introduces  a  compressive  shock  into  the  target  plate.  This 
compressive  shock  traverses  the  target  plate  and  reaches  the  free  surface,  there  producing  an 
unloading  wave  that  travels  back  toward  the  impacted  surface.  The  hypothesis  of  Kanel  et  al. 
was  that  in  K-19  glass,  before  reaching  the  impacted  surface,  the  unloading  wave  encounters  a 
slower  moving  failure  front  (figure  lb).  The  abrupt  change  in  shock  impedance  at  the  failure 
front  causes  a  partial  reflection  of  the  unloading  wave.  The  reflected  wave,  upon  reaching  the 
free  surface,  produces  a  second  plateau  in  the  velocity  signal. 

To  produce  this  hypothesized  phenomenon  in  a  simulation,  a  damage  model  should  contain  two 
features.  First,  the  damage  evolution  equation  should  introduce  one  or  more  time  scales,  thereby 
allowing  the  damage  front  to  lag  the  initial  compressive  shock.  Second,  the  shock  impedance 
should  be  substantially  altered  from  its  pre-damaged  level. 

A  particular  damage  model  that  contains  these  features  was  presented,  initially  in  an  unpublished 
conference  paper  by  Grinfeld  and  Wright  (unpublished),  and  subsequently  in  the  published 
conference  paper  Grinfeld,  Schoenfeld,  and  Wright  (in  press).  This  model  degrades  the  stiffness 
of  an  isotropic  linearly  elastic  material.  The  model  introduces  a  state  variable  that  is  a  measure 
of  damage.  This  state  variable  is  initialized  to  zero  and  then  monotonically  grows  according  to 
an  evolution  equation  that  is  based  on  elastic  strain  energy  and  introduces  a  single  time  scale. 

The  damage  model  also  consists  of  a  degradation  function,  a  function  of  the  damage  state 
variable,  which  is  applied  multipicatively  to  the  material’s  shear  modulus.  The  damage  model  of 
Grinfeld  and  Wright  is  described  in  more  detail  in  section  2. 

The  damage  model  of  Grinfeld  and  Wright  was  installed  into  the  LS-DYNA1  finite  element  (FE) 
code  via  the  user-material  interface  feature  (Livermore  Software  Technology  Corp.,  2003).  This 
implementation  and  the  subsequent  application  of  LS-DYNA  to  problems  of  constant  stretch  rate 
are  described  in  section  3. 

In  section  4,  the  damage  model  is  applied  to  simulations  of  normal  plate-on-plate  impact.  First, 
an  initial  value,  boundary  value  problem  (IYBVP)  is  derived.  Scaling  of  the  IVB  VP  reveals  the 
presence  of  a  dimensionless  group  n,  which  is  a  measure  of  the  material’s  damage  sensitivity. 
The  LS-DYNA  implementation  was  used  to  obtain  solutions  corresponding  to  a  range  of  II 
values.  For  certain  II  values,  the  results  for  the  target’s  free-surface  velocity  showed  a  gradual 
increase  over  time  following  the  initial  arrival  of  the  compressive  shock.  For  other  II  values,  the 

1  LS-DYNA,  which  is  not  an  acronym,  is  a  trademark  of  Livermore  Software  Technology  Corp. 


1 


free-surface  velocity  gradually  decreased  over  time  following  the  initial  arrival  of  the 
compressive  shock.  These  observations  are  summarized,  and  suggested  modifications  of  the 
model  are  presented  in  section  5. 


Figure  1.  The  empirical  evidence  for  failure  waves:  (a)  VISAR  signal  on  the  rear 
surface  of  a  K-19  glass  target  plate  (reproduced  from  Kanel  et  al.,  1992); 
(b)  the  hypothesized  x-t  diagram  for  the  target  plate  (reproduced  from  Brar 
and  Bless,  1992). 
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2.  The  Damage  Model 


The  model  for  damage  in  brittle  materials  was  developed  by  Grinfeld  and  Wright  and  presented 
in  Grinfeld  and  Wright  (unpublished)  and  Grinfeld  et  al.  (in  press).  The  components  of  the 
model  are  described  in  this  section. 


2.1  Pre-damaged  Material 

The  model  is  applied  to  a  material  that  in  its  pre-damaged  state  is  isotropic  and  linearly  elastic. 
Such  a  pre-damaged  material  is  characterized  by  an  elastic  shear  modulus,  //,  Poisson  ratio,  v, 
and  density,  p.  The  material’s  strain  energy  per  unit  volume,  Wo,  is  related  to  the  infinitesimal 
strain  tensor,  e,  by 

^(C)  =  /Ife^ +eySy]  (1) 


in  which  the  summation  convention  on  repeated  indices  applies.  In  Cartesian  coordinates,  the 
components  of  e  are  related  to  those  of  the  displacement  vector  u  and  position  vector  x  by 


1  f  8ut  8uj ' 

2  [a*.  +  dxt  y 


(2) 


The  conditions  of  small  displacement  and  small  strain  have  been  assumed.  Components  of  the 
Cauchy  stress  tensor  <r  are  related  to  Wo  by 


2.2  Damage  Evolution 

The  model  of  Grinfeld  and  Wright  introduces  damage  effects  by  means  of  an  internal  state 
variable  D  that  is  a  function  of  position  and  time.  This  variable  is  initialized  to  zero  throughout 
the  material.  Thereafter,  D  monotonically  increases  according  to  the  evolution  equation 


D 


r^L 

8D 


(4) 


Here,  ///  is  the  Helmholtz  free  energy  per  unit  volume  and  C  is  a  material  constant  with 
dimensions  of  time  x  distance/mass. 
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2.3  Degradation  Function 

Degradation  function  </>(D)  is  used  to  decrease  the  material’s  stiffness  as  damage  accrues.  <p(D) 
relates  the  Helmholtz  free  energy  density  of  the  damaged  material  to  the  elastic  strain  energy 
density  by 

^(e,  D)  =  </>(D)  W0(e)  (5) 


Grinfeld  and  Wright  assumed  (j)  to  be  linearly  related  to  D  according  to 

«t(D)  =  l-(l-«Smi„)D  (6) 

where  <f)min  is  a  second  material  constant  introduced  by  the  model  (C  in  equation  4  was  the  first). 
Equation  6  is  sketched  in  figure  2.  This  dimensionless  constant  tf>min  is  used  to  specify  the 
residual  stiffness  of  fully  damaged  material;  tf>min  satisfies  the  restriction 

o  s  A,„  £  1  (7) 


The  lower  limit  corresponds  to  the  case  of  zero  residual  stiffness  and  the  upper  limit  to  the  case 
of  no  degradation  of  stiffness.  The  combination  of  equations  4,  5,  and  6  yields 


- e  e ..  +  e  .e- 

l-2v  "  JJ  9  v 


(8) 


Coleman  and  Gurtin  (1967)  applied  the  Clausius-Duhem  Inequality  (the  statement  that  the  rate  of 
entropy  production  must  be  non-negative)  to  material  models  involving  one  or  more  internal 
state  variables.  For  the  special  case  of  effectively  isothermal  conditions  (negligible  spatial 
gradients  of  temperature  and  time  rates  of  change  of  temperature)  and  a  single  internal  state 
variable  D,  the  Clausis-Duhem  Inequality 


Equations  4  and  9  lead  to  the  requirement 

C>0  (10) 


requires  that 

dy/ 


dD 

that 


-D<  0 


(9) 
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Figure  2.  Degradation  function  </>(D) . 


2.4  Calculation  of  Stresses 

Coleman  and  Gurtin  (1967)  also  established  a  relationship  between  stress  and  the  Helmholtz  free 
energy  that  follows  as  a  consequence  of  applying  the  Clausius-Duhem  Inequality  to  a  material 
model  involving  internal  state  variables.  In  our  special  case  of  a  single  internal  state  variable  D 
and  small  deformations,  Coleman  and  Gurtin  showed  equation  3  to  be  generalizable  to 


From  equation  5, 


a,i  de„ 


dW 


The  combination  of  equations  1  and  12  yields  the  stress-strain  relation 


cr  =  2<f>(D)/u 


l-2v 


+  G 


in  which  <5y  is  the  Kronecker  delta  function.  The  components  of  the  Cauchy  stress  tensor  are 
then 


^ _ V^xx  +  V(y6yy  +  gzz ^ 


(14a) 


cr  = 

yy 


(14b) 
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crzz  =(j){D) 


l_2vb  V)ezz+Vi.exX+eyy)] 

(14c) 

(14d) 

(Tyz  =</>(D)  2/ieyz 

(14e) 

=<KD)  2//ezx 

(14f> 

This  completes  specification  of  the  model  presented  by  Grinfeld  and  Wright  (unpublished). 


3.  Implementation  Into  LS-DYNA 


3.1  Generalization  to  Large  Strain  and  Large  Displacement 

Prior  to  implementation  of  the  damage  model  into  FE  software,  it  was  first  generalized  to  apply 
to  problems  involving  large  displacements  and  large  strains.  (In  a  proper  application  to  a  brittle 
material,  the  model  is  unlikely  to  encounter  large  strains,  but  large  displacements  and  rotations 
cannot  be  ruled  out.)  The  relevant  continuum  mechanics  concepts  are  explained  further  in 
Malvern  (1969). 


3.1.1  Pre-damaged  Material 


In  the  expression  for  the  elastic  strain  energy  density  function,  the  infinitesimal  strain  tensor  is 
replaced  with  the  Green  strain  tensor  E,  i.e., 


W{](V)  =  M 


- EE  +E..E.. 

\-2v  "  JJ  l}  ,J 


(15) 


where  E  is  defined  in  terms  of  the  deformation  gradient  tensor  F  and  the  identity  tensor  I  by 


E  =  — (Fr-F  -i) 

2V  ’ 


(16) 


In  Cartesian  coordinates,  equation  16  is  expressed  in  terms  of  the  displacement  vector  u  and  the 
material  coordinate  vector  X  by 


**  =  2 


du[ 


Suj  duk  duk 


dX; 


8Xt  dXj , 


(17) 


The  second  Piola-Kirchhoff  stress  S  is  then  computed  from 
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(18) 


S,,=^ 

J  dE, 

and  the  Cauchy  stress  o  is  computed  from 

F  S  Fr 

<T  = - 

detF 


(19) 


3.1.2  Damage  Evolution  and  the  Degradation  Function 

Equations  4  and  6  still  apply.  Helmholtz  free  energy  density  is  now  a  function  of  D  and  E,  i.e., 

V(E,D)  =  <j>(D)W()(E)  (20) 


and  equation  8  is  replaced  with 


— - — EE  +E..E.. 

l-2v  "  JJ  l]  9 


(21) 


3.1.3  Calculation  of  Stresses 

In  order  to  satisfy  the  property  of  frame  indifference,  it  is  convenient  to  work  with  the  second 
Piola-Kirchhoff  stress  S  rather  than  directly  with  o.  Equation  1 1  is  replaced  with 


s„=^- 

J  dE:: 


(22) 


which,  with  equations  1 5  and  20,  yields  the  stress-strain  relation  in  terms  of  the  second  Piola- 
Kirchhoff  stress  and  the  Green  strain. 


S„=2<i>(D)E 


l-2v 


Eh  5  +  E 

kk  ij  ij 


(23) 


The  six  individual  components  of  S  are 

=AD)  [(1  -  v)E  „  +  v(Eyy  +  E„ )]  (24a) 

s„  =^)'fr|:  h  -  V)E„  +  V(£„  +  £„ )]  (24b) 

S„  ^L.[(\-v)E=  +v{Ea+Ey))\  (24c) 

Sxr=^D)-2fiExr  (24d) 

S„=ftD)-2f,E„  (24e) 
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szx  =^pyinEa 

As  with  the  pre-damaged  material,  equation  19  again  determines  the  Cauchy  stress  tensor  a. 


(24f) 


3.2  Implementation  With  the  User  Material  Feature  of  LS-DYNA 

LS-DYNA  offers  a  user  material  interface  whereby  a  user-provided  FORTRAN  subroutine  can 
be  linked  to  the  remainder  of  the  code.  The  user  sets  up  an  array  of  material  constants  (in  this 
case  p,  p,  v,  C,  and  <j>mw)  and  identifies  one  or  more  “history  variables”  (in  this  case  D).  The 
subroutine  is  called  for  every  element  at  every  time  step.  The  code  passes  to  the  subroutine  the 
material  constants,  F  for  the  current  time  step,  and  D  from  the  previous  time  step,  and  the 
subroutine  returns  the  Cauchy  stress  tensor  a. 


The  flow  of  the  calculations  is  as  follows.  Equation  16  is  used  to  compute  E.  Equation  21  is 
integrated  to  update  D  as 


l-2v 


(e)  17(e) 


-EWE 


j] 


i  p(e)  p(e) 


At 


(25) 


The  “(e)”  superscript  indicates  application  to  a  particular  finite  element.  Equation  6  updates  (/), 
equations  (24a-f)  update  S,  and  equation  19  updates  o. 


3.3  Constant  Stretch  Rate  Tension  and  Compression 

In  order  to  study  features  of  the  damage  model  prior  to  its  application  to  normal  plate-on-plate 
impact,  the  model  was  first  applied  to  problems  involving  a  prescribed  constant  stretch  rate. 


3.3.1  Uniaxial  Tension  at  a  Constant  Stretch  Rate 

Figure  3  shows  a  single  eight-node  hexagonal  “brick”  element.  Note  the  material  and  spatial 
coordinate  systems  defined  in  the  figure.  The  “shape  function”  of  this  eight-node  element 
restricts  its  internal  velocity  field  to  the  following  quadratic  form  in  x,  y,  and  z: 

v*  (x,  y,  z,  t)  =  ax  (i t )  +  bx  (t)x  +  cx  (t)y  +  dx  ( t)z  +  ex  ( t)xy  +  /,  (t)xz  +  gx  ( t)yz  (26a) 

vy  ( x ,  y,  z,  t)  =  ay  (t)  +  by  (t)x  +  cy  ( t)y  +  dy  ( t)z  +  ey  ( t)xy  +  fy  ( t)xz  +  gy  ( t)yz  (26b) 

vz  (x,  y,  z,  t )  =  az  {t)  +  bz  (t)x  +  cz  ( t)y  +  dz  ( t)z  +  ez  (t)xy  +  fz  (t)xz  +  gz  ( t)yz  (26c) 

in  which  vx,  vy,  and  vz  are  the  three  components  of  velocity.  The  coefficients  a x  (t),  bx  (t),  ...,g_(t) 
are  functions  only  of  time. 
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At  time  t  =  0,  the  element  in  figure  3  is  a  cube  with  edge  length  Zq.  Thereafter,  a  time-independent 
x-velocity  vo  is  imposed  on  the  four  nodes  of  the  face  X=  Lo.  At  these  same  four  nodes,  the  y-  and 
z-velocities  are  held  at  zero.  At  the  four  nodes  of  the  opposite  face,  defined  by  X=  0,  all  three 
velocity  components  are  held  at  zero.  The  boundary  conditions  of  the  problem  are  therefore 

vx(L0,Y,Z,t)  =  v0t  (27a) 

v,  (0, 7,  Z,  0  =  (0, 7,  Z,  0  =  vz  (0, 7,  Z,  t)  =  vy  (Z0  ,Y,Z,t)  =  vz(L0,Y,Z,t)  =  0  (27b) 

V  7,  Z  e  [0,  Z0  ] 


Once  conditions  27a  and  27b  are  imposed  on  the  shape  functions  of  equations  26a-c,  we  find 
that  motion  throughout  the  element  is  described  by  the  mapping 


x(X,t)  =  X  +  ^-X  = 

Lr, 


y  =  Y 
z  =  Z 


1  +  ^ 

v  A) 


X 


(28a) 

(28b) 

(28c) 


9 


Equations  28a-c  are  exact,  regardless  of  how  fast  one  pulls  the  element.  There  are  no  inertia 
effects  in  this  single-element  problem,  which  makes  it  a  useful  vehicle  for  isolating  effects  of  the 
constitutive  model. 


Thus,  the  deformation  gradient  tensor  is 


F(t)=*  = 
dX 


1  +  ^ 
K 
o 

0 


0  0 

1  0 
0  1 


(29) 


and  the  Green  strain  tensor  is 


E(t)  =  |(FT  -F  -  l)  = 


V 


1+  Vo' 


2  L 


o  J 


0 

0 


0  0 

0  0 
0  0 


(30) 


Note  that  the  strain  is  spatially  constant  throughout  the  element.  If  the  strain  is  small,  then  there 
is  little  difference  between  strain  and  stretch,  i.e., 


Exx(t)  = 


V 


vj 


1  +  - 

V  2L0  j 


vj 


(31) 


and  our  constant  stretch  rate  corresponds  approximately  to  a  constant  strain  rate;  call  it  E . 

p  ~  =  p 

^  XX  —  J  —  ^ 

The  elastic  strain  energy  density  throughout  the  element  is  then 


(32) 


-  EE  +E..E.. 

l-2v  "  ]]  11  lJ 


\ 


PE^{Exx(t)yJAz^.(Ety 

l-2v  l-2v  V  ’ 


(33) 


The  damage  evolution  equation  becomes 

^  =  (1  -  A,  Jc '  w,  (0  =  (1  -  )c .  •  (Et  f 

at  1  -  2v 

The  solution,  assuming  no  initial  damage  at  time  zero,  is  the  cubic  in  time 

D(t)  =  (l  -  </>min  )C  •  •  E2t3 

\  3(i_2v) 


(34) 


(35) 
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However,  this  solution  must  be  modified  to  take  account  of  the  upper  bound  of  1 ,  reached  at 
some  time  t*  determined  from 


£>('*)  =  (!-  Ajc  •  ■  E'-  («*)3  =  1 

3(1  -2v) 


to  be 


t*  = 


1  3(1  -2v) 


n>/3 


(Wmin)C  (1  -  V)// 


-2/3 


(36) 


This  provides  a  time  scale  associated  with  damage  evolution.  The  time  scale  decreases  with 
increasing  strain  rate.  The  two-part  solution  for  D  is  then 


D(t)  = 

The  degradation  function  is 


M  J 


;  0  <t<V 


t>t* 


K, 


;  0  <  t  <  t* 


In  order  to  calculate  stress,  we  return  to  equation  33. 


Then 


From  equation  3 1 , 


so  that 


sxx(t)  = 


w  =£ z^E.e  2 

0  1  -  2v  x* 


Sn=AD)-^r-  =  *D)- 2(1  ^ 

1  -  2v 


s„(o=m-2(!  ?M-Et 

l-2v 

l-(l-^min) 

^min 


2(1 -V)// 
1  -  2v 


v^y 


;  0  <  t  < 


f  >f* 


(37) 


(38) 


(39) 


(40) 
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In  figure  4,  equations  36,  37,  and  40  have  been  evaluated  with  the  use  of  the  material  properties 
in  table  1 .  In  addition,  the  same  problem  was  modeled  with  LS-DYNA,  and  the  numerical 
results  were  virtually  indistinguishable  from  the  analytical  results  plotted  in  figure  4. 


Figure  4.  (a)  The  element  in  figure  3  is  stretched  at  a  constant  positive  (tensile)  rate,  and  (b)  D 
and  Sxx  evolve. 


Table  1.  Parameter  values  used  in  figure  4. 


p  (kg/m3) 

3215. 

n  (GPa) 

193.0 

V 

0.1606 

C  (m-s/kg) 

1.0/1  O'" 

^min 

0.1 

Lo  (m) 

0.01 

v0  (m/s) 

8.0*  I04 

In  figure  5,  equation  40  for  Sxx  is  evaluated  for  the  case  of  no  damage  (C  =  0),  and  the  result  is 
compared  with  the  previous  evaluation  for  C  =  l.Ox  10'5  m-s/kg.  We  see  that  the  damage  model 
has  substantially  reduced  the  stress  levels  attained. 


0.00  0.01  0.02  0.03  0.04  0.05  0.06  0.07  0.08  0.09  0.10 


_ Ms) _ 

Figure  5.  The  solution  to  the  problem  in  figure  3  evaluated  with  and  without  damage. 

3.3.2  Uniaxial  Compression  at  a  Constant  Stretch  Rate 

Next,  the  direction  of  vo  in  figure  3  is  reversed  so  that  a  constant  stretch-rate  uniaxial  compression 
test  is  performed.  Figure  6  shows  the  results  for  D  and  Sxx  vs.  t. 

The  D  results  in  figure  6  are  identical  to  those  in  figure  4  from  the  tension  test.  Equation  35  is 
quadratic  in  strain  rate  and  so  does  not  distinguish  tension  from  compression.  Ultimately,  this 
feature  stems  from  the  hypothesis  of  the  damage  model  that  the  rate  of  damage  growth  is 
proportional  to  the  strain  energy,  which  is  quadratic  in  the  strain  components. 
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Nevertheless,  the  stresses  that  develop  in  the  compression  test  are  properly  compressive.  Note 
that,  for  a  given  absolute  value  of  Exx,  the  value  of  Sxx  in  figure  6  is  equal  in  magnitude  and 
opposite  in  sign  to  its  counterpart  in  figure  4. 


Figure  6.  (a)  The  element  in  figure  3  is  compressed  at  a  constant  negative  stretch  rate,  and  (b)  D 
and  Sxx  evolve. 
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4.  Application  to  Normal  Plate-on-Plate  Impact 


4.1  The  Initial  Value,  Boundary  Value  Problem 

In  normal  plate-on-plate  impact,  a  flyer  plate  is  launched  at  a  certain  velocity  into  a  stationary 
target  plate.  If  the  flyer  and  target  plates  are  composed  of  the  same  material,  the  impact  is  said  to 
be  “symmetric.” 

At  time  t  =  0,  let  the  origin  of  our  stationary  coordinate  system  lie  at  the  center  of  the  impacted 
face  of  the  target  plate.  Let  the  x  axis  be  orthogonal  to  that  surface  and  be  directed  into  the 
target.  That  is,  at  t  =  0,  x  =  0  defines  the  impacted  face  of  the  target  plate  and  x  =  L,  where  L  is 
the  initial  thickness  of  the  target  plate,  defines  the  free  surface,  the  motion  of  which  can  be 
measured  with  the  VIS  AR  technique.  The  y  and  z  axes  of  the  Cartesian  coordinate  system  lie  in 
the  impacted  face  of  the  target  plate. 

In  the  central  portion  of  the  target  plate,  i.e.,  the  region  of  sufficiently  small  y  and  z  so  as  not  yet 
to  be  affected  by  unloading  waves  from  the  lateral  edges,  a  condition  of  uniaxial  strain  is  closely 
satisfied.  Let  ux,  uy,  and  u7  be  the  three  components  of  displacement.  Throughout  the  central 
portion  of  the  target  plate, 


u„  =  u,  =  0 


(41a) 


dux  _  dux 
dy  dz 


(41b) 


The  components  of  the  infinitesimal  strain  tensor  e  are  related  to  displacement  components  by 
equation  2.  The  condition  of  uniaxial  strain  requires  that 


e  =  e  =  e  =  e  = e  =0 

yy  zz  xy  yz  zx 


(42a) 


The  only  component  of  strain  not  identically  zero  is  related  to  the  displacement  field  by 


(42b) 


Equations  14a-f  determine  the  components  of  stress. 


(43a) 


=^zz  =</>(&)■ 


l-2v 


(43b) 


°xy  =  Vyz  =  =  0 


(43c) 
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Equation  8  for  damage  evolution  reduces  to 


^  y  V  t  mm/  i  xx 

8t  l-2v 


To  this  we  add  the  Cauchy  x-momentum  equation 


P 


(44) 


(45) 


The  boundary  conditions  are  sketched  in  figure  7a.  On  the  impacted  face  at  x  =  0,  a  time-inde¬ 
pendent  velocity  vo  is  prescribed.  On  the  free  surface  at  x  =  L,  zero  normal  stress  is  prescribed. 
Homogeneous  initial  conditions  are  imposed  on  displacement,  velocity,  and  damage  D. 


Figure  7.  The  boundary  conditions  imposed  on  the  target  in  the  uniaxial  strain  model  of  normal 
plate-on-plate  impact  (a)  before  and  (b)  after  scaling. 


Equations  42b,  43a,  44,  and  45  are  combined  to  yield  two  coupled  nonlinear  partial  differential 
equations  in  two  unknowns,  ux  and  D.  These  equations  and  the  boundary  and  initial  conditions 
constitute  the  following  IVBYP.  In  all  that  follows,  “DE,”  “BC,”  and  “IC”  denote  “differential 
equations,”  “boundary  conditions,”  and  “initial  conditions.” 


DE 


1  d2ux  d2ux 

'7J~dtr~^2' 


dD  dut  ^ 
dx  dx  , 


(46) 


8D 

dt 


(i-Ajc- 


(i  -  r)/<  r  a«,  y 

l-2v  ^  dx  y 


(47) 
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BC 


IC 


3^ 


(o,0 


K0  ? 


<3x 


M  =  o 


(48) 


ux  (x,0)  =  0  , 


^(x,0)=0,  D(x,0)=0 


(49) 


The  domain  is 

x  e  [0,T]  ,  t  >  0 

The  parameter  co  appearing  in  equation  46  is  defined  by 


co  = 


2(1  -v)fJ. 
(1  -  2  v)p 


(50) 


(51) 


In  an  undamaged,  isotropic,  linearly  elastic  material,  co  would  correspond  to  the  “longitudinal 
wave  speed”  or  “dilatational  wave  speed”  (see  Kolsky  1963,  pp.  10-15). 


4.2  Scaling  the  IVBVP 

The  dimensionless  quantities  x ,  t ,  ux ,  vx,  and  &xx  are  defined  by 


L 


t 


c_o£ 

L  ’ 


v0L  ’ 


<r  = 


PC0v0 


These  are  substituted  into  equations  46  through  50  to  yield  the  scaled  IVBVP 
DE 


r\2  *  r\ 2  * 

o  uv  o  uv 


dt1  3x 


£2 


^ d2ux  dD  duy 
D - ^-  +  - 


a  a  9  /v  /"\  /y 

\ ^  C/X  £/X  (yX  y 


dD 

~di 


=  n  • 


dx 


BC 


d_K 

dt 


(o,f)=l  , 


dx 


(U)=o 


IC 


.  (x,0)  =  0  , 


d_K 

dt 


(x,0)=  0  ,  D(x,0)  =  0 


x  e  [0,1] ,  t  >  0 


(52) 

(53) 

(54) 

(55) 

(56) 

(57) 


in  which 
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(58) 


n  =  (i-(!>,ln)c.^-.i.va2 

2c0 

This  scaled  IVBVP  involves  only  two  dimensionless  parameters,  <f)mm  and  n.  When  we  obtain 
solutions  to  this  IVBVP,  parameters  C,  p,  v,  L,  and  vo  do  not  have  to  be  varied  independently  but 
only  in  the  combination  IT;  II  is  proportional  to  C  and  is  a  dimensionless  measure  of  the 
material’s  damage  sensitivity. 


4.3  Dimensionless  Parameter  II 

For  soda  lime  glass,  Brar  and  Bless  (1992)  give  p  =  2500  kg/m3  and  Co  =  5840  m/s.  The  factor 
( p/2co )  for  this  material  is  0.214  kg-s/m4. 

Choose  (j\ nin  =  0.1.  From  equation  58,  we  have 

n  *  0.193  ^^-C-L-v2 
m 

The  target  plate  thickness  is  typically  a  few  millimeters,  so  L  ~  0.005  m.  One-half  the  flyer  plate 
speed  is  typically  a  few  hundred  meters  per  second,  so  v0  ~  300  m/s.  Thus, 

n  *  86.7  •  C  (59) 

m-s 


4.4  Analytical  Solution  for  the  Case  of  No  Damage  (II  =  0) 

For  n  =  0, 

^  =  0 
8t 


(60) 


and  the  IVBVP  reduces  to 
DE 

d2u0  d2u0 


§(o,?)=i.  §4u')=o 

ot  ox 

IC 

r)u 

M0(x,0)  =  0  ,  — T-(x,0)  =  0 

ot 

A  separation  of  variables  procedure  leads  to  the  Fourier  series  solution 


(61) 


(62) 


(63) 
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8 


u0(xJ)  =  t-—YJ 


1 


7rztt(2n-l)z 


■ sin 


- 1  )n x 

•  C 11/1 

(2n-\  )nt 

1 

<N 

_ l 

1 

<N 

_ 1 

(64) 


This  solution  can  be  obtained  in  a  more  convenient  form  by  means  of  Laplace  transforms.  Let 
l[/(x,/)]  denote  the  Laplace  transform  of  / (x,t) .  This  operation  is  defined  by 


l\f(  x,t  )\  =  |  f  (x,t  )e  st  dt 


The  application  of  this  operation  to  equations  61  and  62  yields  the  transformed  BVP 
DE 

d2U(x,s )  2, 


BC 


U(0,s)  = 


dx2 

1 


-s  U(x,s)  =  0 

8U(l,s) 


=  0 


s'  dx 

in  which  U (x, 5)  denotes  the  Laplace  transform  of  u0 (x, t) . 

The  solution  of  equation  66  subject  to  the  boundary  conditions  67  is 


U(x,  s )  =  \  Lf  “  +  e-s(2-i}  ] - L 

s2i  J  l  +  e-2s 


The  Maclaurin  expansion  of  the  final  factor  is 

1 


1  +  e 


-2s 


=  l-e~2s  +e~4s  +... 


so  that  equation  68  can  be  written  as 

U(x,s)  =  \[e~s*  + 
s 

It  is  easily  verified  that 


e~s( 2-x)  _  e~s(2+x ) 


l[(F  -  Cj  -  c2x )  •  H(f  -  Cj  -  c2x)]  = 


-s(c,+c2x) 


(65) 


(66) 


(67) 


(68) 


(69) 


(70) 


in  which  c\  and  c 2  are  constants  and  H(f  -c,  -  c2x) ,  the  Heaviside  unit  step  function  of  the 
argument,  is  defined  by 


1  ;  Cj  +  c2x  <  t 
0  ;  Cj  +  c2x  >  t 


H(?  -  Cj  -c2x)  = 

Therefore,  the  inverse  Laplace  transform  of  equation  69  is 


(71) 
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u0  ( x,t )  =  (t  -  x)  •  H(f  -  x)  +  (t  +  x  -  2)  •  H(f  +  x  -  2)  -  (t  -  x  -  2)  •  H(f  -  x  -  2)  (72) 


in  which  t  e  [0,3)  • 


Each  term  in  this  d’Alembert  form  of  the  solution  corresponds  to  a  reflection  of  the  shock.  The 
solution  is  truncated  after  three  terms  because  the  phenomenon  depicted  in  figure  1  should  likely 
occur  during  t  e  [1,2]  and  almost  definitely  before  t  =  3 .  In  dimensional  form,  equation  72 
becomes 


u0(x,t )  =  —  {  (c0t  - x)  •  H(c0t  -  x)  +  (c0t  +  x  -  2L)  ■  H(c0t  +  x  -  2L) 
co 

-  ( c0t  -  x  -  2 L)  •  H (c0t  -  x  -  2 L)  } 


(73) 


in  which  t  e  [0,  3L  /  c0 )  . 

4.5  Perturbation  Procedure  for  the  Case  of  Small  II 

If  the  dimensionless  group  II  is  small,  then  an  approximate  solution  to  the  IYBVP  of  equations 
53  through  57  can  be  obtained  if  we  perturb  about  the  equation  72  solution  for  zero  II.  That  is, 
let 

ux  =  Uq  +  ur  n  +  o(n2)  (74a) 


D  = 


Dr  n  +  o(n2) 


(74b) 


If  we  substitute  these  expansions  into  equations  53  through  56  and  collect  terms  of  O(II), 
DE 


r\2  *  rs 2  * 

o  ux  o  ux 


3r 


ax2 


[d 


dx2 


1 

8x  dx 


BC 


IC 


dux 

_  /V 

dt 


q(x,0)  =  0  , 


dDx 

~df 


(o,f)=  0  , 


dux 

~df 


^ du0 

V  5x  J 


dux 

dx 


(u')=o 


(.1,0)  =  0  , 


Dx  (x,0)  =  0 


(75) 

(76) 

(77) 

(78) 


Equations  72,  76,  and  78c  can  be  used  to  determine  D\.  From  equation  72, 

du0(x,t ) 


dx 


=  -H(t-x)+  H(t  +  x - 2)  +  H(t  - x - 2) 


so  that 
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y 

=  H(f  -  x)  -  H(f  +  x  -  2)  -  3  H(f  -  x  -  2) 

J 

The  result  is 

Dx  (x,  t )  =  (t  -  x)  •  H(t  -  x)  -  (t  +  x  -  2)  •  H(f  +  x  -  2)  -  3(f  -  x  -  2)  •  H(f  -  x  -  2)  (79) 

in  which  f  e  [0,3)  •  Our  solution  for  D,  valid  to  O(II),  is  then 

D(x,t)  =  [(£  -  x)  ■  H(f  -  Jc)  -  (f  +  x-2)-H(f  +  x-2)  -  3(f  -  Jc  -  2)  •  H(f  -  x  -  2)]  •  II  +  o(n2 )  (80) 

The  dimensional  form  of  equation  80  is 

2 

D(x,t)  =  (1 —  ^min )  C  •  ~~~  [  (c0t-x)-H(c0t-x)  +  (c0t  +  x-2L)-H(c0t  +  x-2L) 

2c0 

-  3(c0t  -  x  -  2 L)  ■  H (c0t  -  x-  2L)  J 

in  which  t  e  [0,3L/co).  Equation  80  is  evaluated  in  figure  8  for  II  of  0.1  and  0.3.  The  O(II) 
solution  is  linear  in  II,  so  the  solution  in  figure  8b  is  just  a  multiple  of  that  in  8a. 

A  combination  of  equations  72,  75,  and  79  leads  to  an  inhomogeneous  wave  equation  with  a 
known  forcing  function.  The  equation  can  be  solved  for  w, ,  subject  to  the  homogeneous 

boundary  and  initial  conditions  in  equations  77  and  78. 

4.6  LS-DYNA  Results  for  Various  Values  of  II 

The  IVBVP  in  equations  53  through  57  can  be  solved  for  arbitary  II  by  means  of  the  LS-DYNA 
implementation  described  in  section  3.  Figure  9  shows  the  FE  mesh  that  was  employed.  The 
length  between  x  =  0  and  x  =  1  is  divided  into  16,000  eight-node  brick  elements.  The  uniaxial 
displacement  condition  of  equation  41a  was  imposed  by  means  of  constraints  on  the  boundary 
nodes.  The  IVBVP  could  have  been  modeled  with  a  single  row  of  elements;  the  64  rows  shown 
in  figure  9  were  used  to  obtain  good  visibility  in  the  mesh  plots  while  maintaining  an  element 
aspect  ratio  of  one.  Gaussian  quadrature  was  performed  with  one  integration  point  per  element. 
Default  artificial  viscosity  coefficients  were  used:  quadratic  viscosity  Q1  =  1.5,  linear  viscosity 
Q2  =  0.06,  hourglass  coefficient  QH  =  0.1. 


r  du0(x,  t ) 
v  dx 
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Figure  9.  FE  mesh  for  the  target  plate:  (a)  entire  thickness  of  the  plate,  (b)  the  x  =  0  end. 


The  solution  was  obtained  for  II  =  0,  0.1,  0.2,  0.3,  and  0.4,  and  for  (f)mm  =  0.1.  Figure  10  shows 
the  results  for  vx  (1,  t) ,  the  free-surface  normal  velocity  as  a  function  of  time.  The  compressive 
shock  arrives  at  t  =  1 .  For  II  =  0,  the  scaled  velocity  jumps  to  2  and  remains  at  that  level  until 
t  =  3 .  As  II  is  increased  successively  to  0.1,  0.2,  0.3,  and  0.4,  the  free-surface  velocity  jumps  at 
t  =  1  to  ever-decreasing  levels.  Furthermore,  for  II  >  0,  the  free-surface  velocity  does  not 
remain  constant  throughout  the  duration  t  e  (1,3) .  Table  2  compares  vx  (1,  t  =  1.125)  and 
vx(l,t  =  2.000) .  We  see  that 

vx  (1,  t  =  2.000)  >  vx  (1,  i  =  1 . 125) ;  II  =  0. 1,  0.2  (82a) 

vx(l,i  =  2.000)  <vx(l,i  =  1.125);  II  =  0.3,  0.4  (82b) 
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Figure  10.  vx  at  x  =  1  (the  free-surface  normal  velocity)  as  a  function  of  t  and  for 
various  Ft  (<pm\n  =  0.1). 


Table  2.  The  scaled  free-surface  velocity,  vx  at  x  =  1 ,  for  t  -  1.125 ,  t  -  2.000  ,  and  the 
difference  between  the  two  for  various  II  (^min  =  0.1). 


n 

£ 

ii 

ii 

to 

vx(x  =  \,  t  =  2.000) 

vx(l,  2.000)  -vx(l,  1.125) 

0 

2.00000 

2.00000 

0.00000 

0.05 

1.95846 

1.97471 

0.01626 

0.1 

1.91831 

1.94310 

0.02479 

0.11 

1.91051 

1.93596 

0.02545 

0.12 

1.90271 

1.92849 

0.02578 

0.121 

1.90189 

1.92775 

0.02586 

0.122 

1.90115 

1.92701 

0.02586 

0.123 

1.90033 

1.92619 

0.02586 

0.124 

1.89959 

1.92545 

0.02586 

0.125 

1.89885 

1.92471 

0.02586 

0.126 

1.89803 

1.92389 

0.02586 

0.127 

1.89729 

1.92315 

0.02586 

0.128 

1.89647 

1.92233 

0.02586 

0.129 

1.89573 

1.92159 

0.02586 

0.130 

1.89491 

1.92077 

0.02586 

0.131 

1.89417 

1.92003 

0.02586 

0.132 

1.89343 

1.91921 

0.02578 

0.133 

1.89261 

1.91839 

0.02578 

0.14 

1.88727 

1.91273 

0.02545 

0.15 

1.87964 

1.90435 

0.02471 

0.2 

1.84236 

1.85764 

0.01527 

0.3 

1.77184 

1.73859 

-0.03325 

0.4 

1.70665 

1.58810 

-0.11856 
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The  slight  increase  in  velocity  with  time  for  II  =  0.1  and  0.2  indicates  the  arrival  of  a  partially 
reflected  unloading  wave.  However,  this  velocity  rise  is  gradual  over  time  and  the  computed 
free-surface  velocity  does  not  exhibit  the  distinct  second  plateau  of  figure  la. 

Figure  1 1  shows  the  computed  spatial  distributions  of  D  and  ^through  the  target  plate’s  thickness. 

These  profiles  are  plotted  at  time  f  =  1.5,  chosen  to  estimate  the  time  at  which  the  unloading  wave 
should  encounter  the  damage  front  according  to  the  hypothesis  sketched  in  figure  lb.  In  figure  11, 

for  a  given  value  of  n,  the  D  and  ^profiles  display  a  resemblance  because  of  the  linearity  of  their 
relationship  in  equation  6.  Note  that  changes  in  (/)  are  broadly  distributed  across  the  thickness  and 
do  not  occur  abruptly  at  a  specific  location.  As  a  consequence,  increases  in  vx  at  x  =  1  for  n  =  0. 1 
and  0.2  occur  only  gradually  with  time  according  to  this  model. 


<pmin  =  0.1 

'  min 

tc0/L=  1.50 


Figure  11.  D  and  (j>  as  functions  of  x  at  t  =  1.5  and  for  various  fl  (<j)mm  =  0.1 ). 


The  solutions  for  n  =  0.1  and  for  n  =  0.3  are  examined  more  closely  in  figures  A-l  through  A-9 
of  appendix  A  and  figures  B-l  through  B-9  of  appendix  B,  respectively. 

Figures  A-l  through  A-3  and  B-l  through  B-3  show  results  for  scaled  stress  <jxx .  One  effect  of 
damage  has  been  to  decrease  the  stress-wave  speed.  For  both  n  =  0.1  and  0.3,  the  compressive 

loading  wave  during  t  e  (0,1)  traverses  the  undamaged  target  plate  with  a  scaled  speed  very 
nearly  equal  to  1,  but  during  t>  1,  the  unloading  wave  travels  back  into  the  damaged  material 
with  a  speed  smaller  than  1 .  Furthermore,  this  wave  speed  reduction  is  more  pronounced  at  n  = 
0.3  than  at  n  =  0.1.  A  second  effect  of  damage  has  been  to  decrease  the  amplitude  of  the  stress 
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signal,  and  this  effect  is  also  more  pronounced  at  II  =  0.3  than  at  II  =  0.1.  In  figure  B-l,  note 

that  the  blue  color  lightens  in  the  wake  of  the  stress  wave  as  early  as  t  =  0.75 ,  indicating  a 
substantial  decrease  (in  absolute  value)  from  the  level  of-1  that  would  pertain  without  damage. 

Figures  A-l  and  A-2  and  figures  B-l  and  B-2  show  through-thickness  profiles  of  axx  at  scaled¬ 
time  intervals  of  0. 125.  These  figures  further  document  the  erosion  over  time  of  the  absolute 
value  of  the  compressive  stress. 

Figures  A-4  through  A-6  for  FI  =  0.1  and  B-4  through  B-6  for  FI  =  0.3  show  results  for  D.  These 
figures  show  that  at  a  given  x ,  location  D  increases  monotonically  with  time.  Figures  A-5  and 
B-5,  which  contain  through-thickness  profiles  of  D  at  fixed  times,  can  be  compared  with  figures 
8a  and  8b,  respectively:  the  perturbation  solution  truncated  at  O(II)  is  useful  at  II  =  0.1  but 
grossly  underpredicts  D  at  II  =  0.3 .  In  figures  A-5  and  B-5,  we  see  that  the  through-thickness 
profiles  of  D  do  not  exhibit  steep  fronts — a  feature  that  was  also  noted  in  figure  1 1 . 

Nevertheless,  if  we  track  the  location  at  which  D  equals  a  specific  constant  between  0  and  1,  that 
location  drifts  to  ever-larger  x  values  over  time  until  the  unloading  wave  arrives  to  relieve  the 
local  strain  energy.  This  is  clearly  seen  in  the  x-t  plots  in  figures  A-6  and  B-6.  These  figures 
also  display  the  location  of  the  stress-wave  front,  which  is  identified  with  the  location  of  the 
maximum  value  of  \daxx/dx\.  Finally,  figures  A-6  and  B-6  also  depict  the  reduction  in  the 

stress-front’s  speed  for  t  >  1  by  displaying  the  (L  -  c0 1 )  line. 

Figures  A-7  through  A-9  and  B-7  through  B-9  present  results  for  scaled  velocity  vx .  Figures  A-7 
and  B-7  show  through-thickness  profiles  at  t  intervals  of  0.125  during  t  e  (0,1) .  Figures  A-8, 
A-9,  B-8,  and  B.9  show  such  profiles  during  t  e  (1,2) .  Figures  A-9  and  B-9  are  enlargements  of 

A-8  and  B-8,  respectively,  in  the  vicinity  of  the  free  surface  at  x  =  1 .  Figures  A-9  and  B-9  clearly 
show  the  increase  in  vx  at  x  =  1  with  increasing  time  for  FI  =  0.1  and  its  decrease  for  FI  =  0.3 . 

The  IVBVP  in  equations  53  through  57  was  then  solved  for  additional  values  of  FI.  Figure  12 
and  table  2  present  our  results  for  the  dependence  of  [vT  (x  =  1  ,t  =  2.000)  -  vx (x  =  1,  t  =1.1 25)] 

on  II.  This  quantity  attains  its  maximum  value  at  FI  =  0.126 .  The  solution  for  this  value  of  FI  is 
displayed  in  figures  C-l  through  C-7  of  appendix  C.  This  solution  is  qualitatively  similar  to  that 
obtained  for  n  =  0. 1 . 
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5.  Conclusions 


The  damage  model  for  brittle  materials  presented  in  Grinfeld  and  Wright  (unpublished)  and 
subsequently  in  Grinfeld  et  al.  (in  press)  was  installed  in  the  LS-DYNA  FE  code.  This  damage 
model  introduces  a  damage  state  variable  D.  The  model  consists  of  an  evolution  equation  for  D 
and  a  degradation  function  of  D.  The  degradation  function  is  applied  multiplicatively  to  the 
elastic  shear  modulus  of  an  isotropic  linearly  elastic  material;  Poisson  ratio  is  unaltered. 

The  damage  model  was  initially  applied  to  the  constant  stretch-rate  uniaxial  tension  and 
compression  of  a  single  eight-node  brick  finite  element.  We  saw  D  and  axial  stress  increase  as  a 
cubic  and  a  quartic,  respectively,  with  time.  We  also  saw  D  attain  the  same  value  for  a  given 
absolute  value  of  stretch,  regardless  whether  that  stretch  involved  tension  or  compression.  Since 
the  evolution  equation  is  based  on  elastic  strain  energy  (a  quadratic  function  of  the  strain 
components),  the  damage  model  does  not  distinguish  damage  contributions  from  tension  and 
compression.  This  issue  can  be  addressed  in  a  future  refinement. 

The  model  was  then  applied  to  simulations  of  normal  plate-on-plate  impact.  We  scaled  the 
IVBVP  and  found  the  occurrence  of  only  two  dimensionless  parameters,  and  II.  The  latter 
is  a  combination  of  pre-damaged  elastic  properties,  material  properties  introduced  by  the  damage 
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model,  target  geometry,  and  impact  speed.  II  is  a  measure  of  the  material’s  damage  sensitivity  to 
elastic  strain  energy.  The  LS-DYNA  implementation  of  the  Grinfeld-Wright  model  was  used  to 
obtain  solutions  for  the  normal  plate-on-plate  impact  problem  for  II  =  0.1,  0.2,  0.3,  and  0.4.  For 
smaller  values  of  II,  a  procedure  for  perturbing  the  II  =  0  (no  damage)  solution  was  outlined. 

The  LS-DYNA  solutions  for  II  =  0.1,  0.2,  0.3,  and  0.4  were  examined  with  regard  to  the  failure 
wave  hypothesis  sketched  in  figure  1.  For  II  =  0.1  and  0.2,  the  normal  velocity  of  the  target’s 
free  surface  was  seen  to  gradually  increase  slightly  from  its  level  associated  with  the  initial 
arrival  of  the  compressive  shock.  For  II  =  0.3  and  0.4,  the  free-surface  velocity  was  seen  to 
gradually  decrease  from  its  level  associated  with  the  initial  arrival  of  the  compressive  shock. 

The  solutions  for  damage  and  stress  fields  within  the  target  revealed  two  competing  mechanisms. 
The  evolution  of  damage  led  to  stiffness  gradients  within  the  target,  which  did  indeed  produce 
gradual  reflections  of  the  unloading  wave  back  to  the  free  surface.  However,  damage  evolution 
also  led  to  a  reduction  of  the  compressive  stress  field,  which  acted  to  reduce  the  free-surface 
velocity. 

A  proposed  change  in  the  Grinfeld-Wright  model  that  might  lead  to  the  distinct  recompression 
plateau  in  figure  lb  is  to  replace  the  linear  $X> )  relationship  of  equation  6.  A  more  abrupt 
change  of  (f)  with  D  may  achieve  both  goals  of  a  more  abrupt  reflection  of  the  loading  wave  and  a 
smaller  degradation  in  the  stress  amplitude  of  the  unloading  wave  before  the  unloading  wave 
encounters  the  abrupt  change  in  (/). 
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Appendix  A.  LS-DYNA  Results  for  II  =  0.1 
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Figure  A-2.  <7XX  as  functions  of  X  for  various  t  e  [0,1]  and  for  fl  =0.1  .  ( (f>mm  =0.1.) 
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Figure  A-3.  (7XX  as  functions  of  X  for  various  t  e  [1,2]  and  for  fl  =0.1  .  ( (fimm  =0.1.) 
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Figure  A-4.  Contours  of  D  across  the  target  plate  for  II  =  0. 1  and  (j)mm  =0.1. 
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Figure  A-5.  D  as  functions  of  X  for  various  t  and  for  fl  =  0. 1 .  ( <f>mm  =0.1.) 
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Figure  A-7.  vx  as  functions  of  X  for  various  t  e  [0,1]  and  for  II  =  0. 1 .  ( (f>mm  =0.1.) 
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Figure  A-8.  vx  as  functions  of  X  for  various  t  e  [1,2]  and  for  II  =  0. 1 .  ( (f>mm  =0.1.) 


38 


n  =  o.i,  i*u,  =  o.i 


c0t/L  =  1.000 
c0t!L  =  1.125 
c0t/L=  1.250 
c0  f  /  L  =  1 .375 
c0t/L=  1.500 
c0  f  /  L  =  1 .625 
c0t/L=  1.750 
c0  f  /  L  =  1 .875 
c0t/L  =  2.000 


Figure  A-9.  Enlargement  of  vx  as  functions  of  X  for  various  t  e  [1,2]  and  for  El  =  0. 1 .  ( (f>mm  =  0.1 .) 
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Appendix  B.  LS-DYNA  Results  for  II  =  0.3 
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Figure  B-l.  Contours  of  Oxx  across  the  target  plate  for  II  =  0.3  and  <pmm  =  0.1 . 
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Figure  B-2.  <TXX  as  functions  of  X  for  various  t  e  [0,1]  and  for  II  =  0.3 .  ( <pmm  =  0.1 .) 
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Figure  B-3.  (7XX  as  functions  of  A'  for  various  t  e  [1,2]  and  for  II  =  0.3 .  ( (f)mm  =  0.1 .) 
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Figure  B-7.  vx  as  functions  of  X  for  various  t  e  [0,1]  and  for  El  =  0.3 .  ( (f)mm  =0.1.) 
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Figure  B-8.  vx  as  functions  of  X  for  various  £  e  [1,2]  and  for  El  =  0.3 .  ( (f>mm  =0.1.) 
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Figure  B-9.  Enlargement  of  vx  as  functions  of  X  for  various  t  e  [1,2]  and  for  El  =  0.3 .  ( (f>mm  =0.1.) 
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Appendix  C.  LS-DYNA  Results  for  II  =  0.126 
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Figure  C-l.  CJXX  as  functions  of  X  for  various  t  e  [0,1]  and  for  II  =  0. 126  .  ( (f>mm  =0.1.) 
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Figure  C-2.  <7XX  as  functions  of  A'  for  various  t  e  [1,2]  and  for  II  =  0. 126  .  ( (f>mm  =0.1.) 
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Figure  C-5.  vx  as  functions  of  X  for  various  t  e  [0,1]  and  for  II  =  0. 126  .  ( (f)mm  =0.1.) 
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Figure  C-6.  vx  as  functions  of  X  for  various  t  e  [1,2]  and  for  fl  =  0.126  .  ( (/)mm  =  0.1 .) 
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FRANCE 

2  RAFAEL  BALLISTICS  CTR 

M  MAYSELESS  ZVI  ROSENBERG 
PO  BOX  2250 
HAIFA  3 1021 
ISRAEL 

1  UNIV  OF  WATERLOO 
MECHANICAL  ENGNRG 
M  J  WORSWICK 
200  UNIVERSITY  AVE 
WEST  WATERLOO 
ONTARIO  N2L3G1 
CANADA 

1  ARL  ERO 
S  SAMPATH 
AERO  MECH  ENG 
223  MARYLEBONE  RD 
LONDON  NW1  5TH 
UNITED  KINGDOM 

1  AMC  SCI  &  TECH  CTR 
EUROPE 

T  J  MULKERN 
POSTFACH  81 
55247  MAINZ  KASTEL 
GERMANY 

2  FRAUNHOFER-INSTITUT 
ERN  S  T -M  ACH-IN  S  TITUT 

V  HOHLER  E  STRASSBURGER 
ECKERSTRASSE  4 
79104  FREIBURG 
GERMANY 
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